Full-time dynamics of modulational instability in spinor Bose— Einstein condensates 
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We describe the full-time dynamics of modulational instability in F = 1 spinor Bose-Einstein 
condensates for the case of the integrable three-component model associated with the matrix non- 
linear Schrodinger equation. We obtain an exact homoclinic solution of this model by employing 
the dressing method which we generalize to the case of the higher-rank projectors. This homoclinic 
solution describes the development of modulational instability beyond the linear regime, and we 
show that the modulational instability demonstrates the reversal property when the growth of the 
modulated amplitude is changed by its exponential decay. 
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I. INTRODUCTION 

Spinor Bosc-Einstein condensate (BEC) with an opti- 
cal confinement represents a unique macroscopic system 
with the spin degrees of freedom [Hi]. The interplay 
between the mean-field effective nonlinearities of three- 
component matter waves and their spin properties pro- 
duce many interesting phenomena such as the spin mix- 
ing H , as well as the formation of spin domains [3, H[ and 
spin textures [!, [f| . Various properties of the spinor BEC 
have been analyzed theoretically @, H, ■ The ground 
state of the spinor BEC with the hyperfine spin F = 1 
can be either ferromagnetic (maximum spin projection) 
or polar (zero spin projection). It was shown in Ref. [10( 
within the linear stability analysis of the spinor conden- 
sate model that the ferromagnetic phase of the conden- 
sate can experience instability for large enough densities 
of atoms, while the polar phase remains always modula- 
tionally stable. 

Wadati and co-authors demonstrated that the 
three-component nonlinear equations describing the evo- 
lution of the F = 1 BEC can be reduced, under special 
constraints imposed on the condensate parameters, to 
the completely integrable matrix nonlinear Schrodinger 
(NLS) equation [12|]. Both bright and dark three- 
component BEC solitons have been found in the frame- 
work of this model [H [3 El (M E3 ■ 

As regards the linear stability analysis presented in 
Ref. [101 , only an initial (linear) stage of the perturba- 
tion development can be explored by this method which 
predicts the exponential growth of the modulation fre- 
quency sidebands for some conditions, i.e., it describes 
the conditions of modulational instability (MI). A physi- 
cal mechanism behind the MI is the parametric coupling 
between the spin degrees of freedom which leads to a pop- 
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ulation transfer between the spin components. To study 
the long-time evolution of instabilities, numerical meth- 
ods are used as a rule. For the scalar NLS equation, the 
problem of the long-time evolution of the MI was stud- 
ied by the truncation of the original model to a finite 
number of modes (as usually, the three-mode approxi- 
mation) [l8j]. More complete analysis of the long-time 
MI dynamics [HI, H(| is based on a linear constraint im- 
posed on the real and imaginary parts of solutions of the 
scalar NLS equation, and it allows one to find a class of 
three-parameter solutions sharing this property. Among 
the solutions found in such a way, a special solution de- 
scribes the development of MI beyond the linear regime, 
and it is identified as a homoclinic orbit separating two 
qualitatively different types of periodic solutions. A sim- 
ilar result was obtained by means of the Darboux trans- 
formation with the plane wave as a 'seed' solution [2lj |. 
Following terminology of Ref. [13] , the full-time dynam- 
ics represents a nonlocal view of the MI development over 
a long time interval. 

A homoclinic orbit is a trajectory of a dynamical sys- 
tem that tends to the same manifold (fixed point, pe- 
riodic orbit, etc.) as time tends to ±oo. The exis- 
tence of homoclinic solutions serves as an indicator of 
chaotic regimes in a perturbed deterministic system. For 
nonlinear wave systems described by partial differential 
equations, the complete understanding of the homoclinic 
structures in the infinite-dimensional phase space is far 
from being available at present. On the other hand, the 
unique features of the integrable nonlinear wave equa- 
tions admit essentially more deep insight into this prob- 
lem. Extended reviews of analytical and numerical meth- 
ods for obtaining homoclinic orbits for the scalar NLS 
and sine-Gordon equations are given in Refs. [23l. l24j . 

The aim of our paper is twofold. First, we derive a 
homoclinic solution of the matrix NLS equation. Sec- 
ond, using these analytic results, we present the exact 
solution of the problem of the long-time evolution of the 
modulationally unstable F = 1 BEC in the case when it 
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is described by the integrable model. 

To find homoclinic solutions, we do not impose ad hoc 
constraints on the form of solutions. Instead, we use a 
kind of dressing procedure, well known in the soliton the- 
ory [25| , which was proposed recently as a systematic tool 
to generate exact homoclinic solutions of i nteg rable non- 
linear equations with periodic boundaries |26J. A dress- 
ing factor being the main technical ingredient in this ap- 
proach contains a projector which determines the coor- 
dinate dependence of the homoclinic solution. It should 
be pointed out that for all known homoclinic solutions 
obtained up to now for various nonlinear equations (see, 
e.g., Refs. [H, this projector has rank 1. A cru- 

cial feature of the matrix NLS equation consists in the 
fact that the corresponding dressing factor incorporates 
the rank 2 projector. In terms of the soliton theory it 
corresponds to multiple zeros of the scattering matrix 
coefficients (or multiple zeros of the associated Riemann- 
Hilbert problem). Notice that the case of multiple zeros 
cannot be treated as a coalescence of simple zeros [28j . 
Accordingly, we modify the definition of the dressing fac- 
tor for the case of the matrix NLS equation and obtain 
the first example of the matrix homoclinic orbit and, as 
a result, the complete description of the MI evolution in 
the integrable spinor BEC model. 

The paper is organized as follows. In Sec. II we de- 
scribe the integrable F = 1 BEC model. The method 
for obtaining homoclinic solutions for integrable nonlin- 
ear equations valid for higher rank projectors is outlined 
in Sec. III. Section IV is devoted to the explicit deriva- 
tion of the homoclinic solution for the matrix NLS equa- 
tion and presents the main results of our paper. The 
homoclinic solution describes the temporal evolution of 
linearly unstable modes. We show that the MI has a 
reversal property - the initial-wave profile is recovered 
after a sufficiently long time. Hence, the term 'side-band 
instability' refers in fact to only the linear stage of the 
instability development. Section V concludes the paper. 



II. MODEL 

We consider an effective one-dimensional BEC trapped 
in a pencil-shaped region elongated in the x direction 
and tightly confined in the transversal directions. The 
assembly of atoms in the hyperfine spin F — 1 state 
is described by a vector order parameter $(x,t) = 
(<f>+(x, t),$o(x, t), <&-(x, t)) T , where its components cor- 
respond to three values of the spin projection mp — 
1,0,-1. The functions <&± and $n obey a system of 
coupled Gross-Pitaevskii equations |29| 
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where the constant parameters cq = (go + 2g 2 )/3 and 
c 2 = {92 — <7o)/3 control the spin-independent and spin- 
dependent interaction, respectively. The coupling con- 
stant gf (/ = 0, 2) is given in terms of the s-wave scat- 
tering length a/ in the channel with the total hyperfine 
spin /, 
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Here a± is the size of the transverse ground state, m is 
the atom mass, and C = —£(1/2) « 1.46. 

It was noted in [U that Eqs. flO} are reduced to an 
integrable system under the constraint 



co = c 2 



-c < 0. 



(2.2) 



The negative c 2 means that we consider the ferromag- 
netic ground state of the spinor BEC with attractive in- 
teractions. The condition (|2.2[) . being written in terms 
of gf as 2go = — g 2 > 0, imposes a constraint on the 
scattering lengths: a± = 3Caoa 2 /(2ao + a 2 ). Redefin- 
ing the function $ as <l — * (4>+, V%4>0i 4>-) T , normalizing 
the coordinates as t — > (c/H)t and x — > (y/2mc/h)x, and 
accounting for the constraint (12. 2\i , we obtain a reduced 
system of equations in a dimensionless form: 

id t( j> ± +dl(j) ± +2 (|0±| 2 + 2|0 O | 2 ) 0±+2(^ = , (2.3) 

id t ^+d 2 x 4>o+2 (|0+| 2 + \M 2 + I'M 2 ) = . 

After arranging the components <f>± and (f>o into a 2 x 2 
matrix Q = ^ ^ + ^° J i we can easily see that Eqs. 
take the form of the integrable matrix NLS equation 

id t Q + d 2 x Q + 2QQ t Q = . (2.4) 



The matrix NLS equation (|2.4[) possesses the Lax repre- 
sentation with the 4x4 matrices U and V of the form 

m 



U = ikA + 0, A = diag(-l, -1, 1, 1), Q 



Q 
-Q f 

(2.5) 

V = 2ik 2 k + 2kQ + i( Q ^ _gtg) . ( 2 - 6 ) 
k is a spectral parameter. 



III. METHOD 

We are interested in periodic solutions of Eqs. (|2.3[) 
(or (|2.4p ) with a spatial period L, Q(x + L,t) = Q(x,t). 
Hence, the Floquct theory should be applied to analyze 
the spectral problem 
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The fundamental solution M(x, k) of Eq. (|3.ip is fixed 
by the condition M(0, k) = I, I is the unit 4x4 matrix. 
Then we define a transfer matrix T(k) as the fundamen- 
tal solution in the point x — L, T(k) = M(L,k). Di- 
agonalization of the transfer matrix determines a matrix 
R, 

R^T^R = diag (e lmiL , e lm * L ) = A(L, k), 

and produces the Floquet multipliers exp(zm :) L) with the 
Floquct exponents uij , j — 1, . . . , 4. The Floquet spec- 
trum is a set of all k for which the transfer matrix T(k) 
has the eigenvalues on the unit circle. 

The next step is a determination of a Bloch solution 
X of Eq. (|3.ip as x = MR which obeys the property 
x(x + L,k) = x(x, k)A(L, k), specific for the Bloch-type 
solutions. The Bloch eigenfunctions of the periodic spec- 
tral problem (|3.1[) play the role of the Jost solutions of 
the spectral problem with a decreasing potential. 

Among the points of the Floquet spectrum we will dis- 
tinguish the so called double points [2^]. Double points 
arc those values of k for which the Floquet exponents 
rrij differ in multipliers of 2n/L or, in other words, the 
Floquet multipliers are degenerate. We will be especially 
interested in complex double points which indicate lin- 
earized instability of solutions of Eq. (|2.3[) and label or- 
bits homoclinic to hyperbolic fixed points in the phase 
space of a nonlinear system. Note that the term "dou- 
ble" in the context of the Floquet spectrum refers to the 
algebraic multiplicity of a point of the spectrum and has 
no relation to the multiplicity of zeros we have spoken 
about in the Introduction. Real double points are asso- 
ciated with stable modes. 

Suppose we know explicitly a Bloch solution xo of the 
spectral problem xox = UqXo with the matrix Uq (|2.5p 
whose entries contain the known solution Qo of Eq. (|2.4[) . 
Then we dress the solution xo by applying the dressing 
factor D(x,t,k), x = Dxo, and % is a new solution of 
the spectral problem with new matrix U = DUqD^ 1 + 
D X D~ X . The dressing factor has the form 



N h — h* 
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where P s is a projector, P s 2 = P. 
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Here k s , s = 1, . . . ,N are complex double points of the 
Floquet spectrum and r s is the rank of the projector P s . 
The four-component ket- and bra-vectors \n; s) and (I; s\ 
are the column and row arrays, respectively. Hence, |n; s) 
is a four vector related with the sth complex double point 
k a and obtained by applying the Bloch solution xo{k s ) 
to a constant vector \q; s), 



Exactly A s *> vectors \q;s), and hence r^ s ' vectors \n;s), 
correspond to the complex double point k s . The sum- 
mation in Eq. (|3.2p is taken over all N complex double 
points, while that in Eq. (|3.3[) is performed over the 
r s -dimensional space of vectors \n; s) produced in accor- 
dance with Eq. (|3.4|) . Then a new solution of the matrix 
NLS equation is written as 



N 



Q = Qo + X>s-O a > p *]- 



(3.5) 



For the rank 1 projectors these formulas reduce to the 
known ones [26j |. 

Note the essential difference in applications of the 
dressing procedure between the soliton theory and the 
periodic wave theory. Indeed, the parameters k s are free 
in the standard use of the dressing method and, in any 
case, they do note relate with the seed solution Qq. On 
the contrary, our approach demands to choose k s as the 
complex double points of the Floquet spectrum of the 
spectral problem (|3.f [) for the seed solution xo- 

Hence, they are the complex double points k s and the 
projectors P s that completely determine new solution. 
In the next section the above method will be used to 
generate homoclinic solution of the spin 1 BEC model 
(|2.3p and hence to reveal the long-time dynamics of the 
MI in this model. 



IV. RESULTS 

We begin with a spatially homogeneous continuous 
wave solution of Eq. (|2.3j) with components 



if = 



= ibe 



-if.it 



(4.1) 



as the seed solution to be dressed. Here a and b are 
real constant amplitudes which determine a population 
of each spin component, and the chemical potential fi is 
given by /i = —2 (a 2 + b 2 ). Note the fixed ir/2 phase 



difference between the components (j>^ and ■ The 
same phase locking property is an inherent feature of the 
nonintegrable model (|2.1[) as well [l(J. We could start 
with more general representation of plane waves but the 
structure of Eqs. (|2.3[) and the Galilean invariance re- 
duce it to the form (|4.1[) . Then we consider the spectral 
problem (|3.1[) with the matrix Uq containing the plane 
waves (|4.1[) as the potential Q : 



Uo = 



—ik 


ibe iflt 
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ibe~ iflt ae- 1 ^ 
ik 
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(4.2) 



\n;s) = xo(x,t,k s )\q; s). 



(3.4) 



The fundamental solution of the spectral problem with 
the matrix Uq is explicitly found: 
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( cospx + i(k/p) sin pa; (a/p) sinpxe~ vt i(b/p)s\a.pxe~ l ^ t \ 

cos pa; + i(k/p) s'mpx i(b/p) sin pxe~ lfJ,t (a/p)sinpxe~ l ^ t 

— (a/p) sinpxe v * i(b / p) sinpxe 1 ^ 1 cospx — i{k/p) sin px 



\ i(b/p)sinpxe 
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-(a/p) s'mpxe^ 1 







detM = 1, 



ipx — i(k/p) sinpx J 



(4-3) 



where p 2 = a 2 + b 2 + k 2 . Diagonalization of the transfer 
matrix T(k) = M(L, k) is performed by the matrix R 



which has the form 
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where dj are time dependent. As a result, 

R- 1 T(k)R = A(L, fc) = diag(e" lpL , e- jpi , e lpi , e lpL ) . 

Therefore, the Floquet exponents are written as 

mi = — p, m 2 = —p, rri3 = p, = p (4.5) 

and have the multiplicity 2. Then we obtain the seed 
Bloch solution \o = MR in the form 



Xo = exp(i-t) 
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exp(zpAa; + 2ikpAt), 



(4.6) 

where the parameters dj (t) entering Eq. (I4.4[) have been 
determined from the second Lax equation xot — VqXo 
with the matrix Vq (|2.6p depending on the seed continu- 
ous wave (14.11): 



di = ^10 cxp(— 2ik 2 t), d 2 = d 2 oexp(2ik 2 t) 7 



d 3 =d 30 exp 



—fit — 2ikpt 



c?4 = c?4o exp 



—/it + 2ikpt 



Here djo are integration constants. 

Now we proceed to finding the complex double points. 
Following Ref. [13] , we seek for double points as a differ- 
ence between two Floquet exponents mi and m 3 (|4.5[) : 



m 3 = mi + 6 a 



2n 



-s, 



s = ±l,±2,.. 



This gives 

fc s = ±i y / a 2 + b 2 - (irs/L) 2 , 
and 

k s = ±^(n S /L) 2 -(a 2 + b 2 ) } 



if a 2 + b 2 > (irs/L) 2 
(4.7) 



if a 2 + b 2 < (ns/L) 2 . 

(4.8) 

Hence, for given amplitudes a and b and period L the 
double points are arranged into infinite number of real 
double points (|4.8p situated on the real axis in the k 
plane, and a finite number of complex double points 
(|4.7p lying on the imaginary axis within the interval 
(-Wa 2 + b 2 ,iVa 2 + b 2 ). 

Let us choose in the following the amplitudes and pe- 
riod in such a way that to obtain the single complex 
double point ki (and hence — k\). It means \Ja 2 + b 2 > 
(ir/L) but \Ja 2 + b 2 < (2n/L). For this choice the only 
rank 2 projector P = Pj has the form 



P = 



1 



n\l) 



ki - ki ^ 

1 nl=l 



To simplify notations, we write \n) instead of \n; 1). Since 
D is a 2 x 2 matrix, we easily obtain the following ex- 
pression for the projector: 



P=i[(2|2)|l)(l| 



(1|2)|1)(2|-(2|1)|2)(1| 



(1|1>|2)(2|], 
(4.9) 



D = (l|l)<2|2)-(1|2)<2|1), 
where the vectors |1) and |2) are determined as 
|l)=Xo(fei)k), |2)= X o(*i)|r). 
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Here \q) and \r) are linearly independent constant vectors 
with the components qj and rj, j — 1, . . . ,4. Then after 
rather lengthy but straightforward algebraic calculation 
in accordance with Eqs. (|4.9[) and (|3.5[) taken for s = 1, 
we explicitly obtain new solutions of the integrable spin 
1 BEC model (1231). 



ae -v* ( l + 2i— sin^> ) , (4.10) 



„ =- ih, ''" [ l + 2^sinV 



which at the same time represent components of the ma- 
trix homoclinic orbit of the matrix NLS equation (|2.4ll . 
Here 

A = cosh2r — cos2psin 2 ip + 27 cosh r sin p sin ^ + -7' ! , 



(441) 



_B = sinh2rcos^ + i cosh 2r sin -0 — icos2psin-0 



+ i 7 [( M _e T + < e^e 4 ' - (/i* e T + M+e-^e-v 



+ ^7 2 ( sinV> — — cos ib cos 023 I , (442) 
2 V a 



i?o = sinh 2r cos ?p + i cosh 2t sin ^ — i cos 2p sin V> 




FIG. 1: Full-time evolution of the <f>+ (and (/>-) component 
due to modulational instability. The parameters are a = 1, 
b = 2, L = tt/2, q 2 = tt/3, a 3 = vr/ 4 > M = 1- 



i/j- = 1 ± ie 



^sin ip 



The solutions (|4.10p are indeed homoclinic to the plane 
waves (|4.ip . Calculation of the asymptotics of (j>± and 0o 
as t — > ±00 gives 



+ -7 [(^_e r + ^*e- T )e v - «e r + ^e^e^'] 



i 2 ( ■ , a 1 

--7 ^S1H l/> — — COS l/> COS «23 



(4.13) 



The constants djo have been incorporated into the con- 
stant components qj and rj of the vectors \q) and \r). If 
we denote definite combinations of these components as 

ei = qir 2 - q^i, e 2 = qir 3 - q%r\, e 3 = q x r± - q 4 ri, 
C4 = 93'"4-g47'3, e 3 - = |e-,|e Mj , ay; = a 3 -ai, |e 2 | = |e 3 |, 



then the notations used in (14. lip . (|4.12|) and (14. 43)) are 
as follows: 



fci = ik , pi = J a 2 + b 2 + k\ = — , 



,:</; _ Pi + ik 
~ Va 2 + b 2 : 



In other words, these solutions reproduce in the limit 
t — > ±00 the seed plane waves up to a constant phase, 
as should be for the homoclinic orbit. Note that the 
nonlinear MI for the spin 1 condensate but for different 
phases of the seed wave components was studied by the 
Darboux transformation in Ref. [14j. 

Figures [T]and [2] illustrating the solution (|4. 10[) demon- 
strate typical development of the continuous wave pertur- 
bation within three periods in x. We see that the stage 
of the exponential growth of instabilities revealed by the 
linear stability analysis transforms to the exponential de- 
creasing with emergence of localized structures. Hence, 
the full-time evolution of MI for the integrable F = 1 
BEC model demonstrates the reversal property, such as 
the Fermi-Pasta-Ulam process [3(J: the phase trajectory 
of the system returns to the initial one which corresponds 
to the continuous waves (|4.1[) . For chosen parameters the 
growth and decrease development of the component 4>o 
is more pronounced than that of <j>± . 



t = Akopit + 1 , p = 2p 1 x-a 13 , a 13 = a 34 , 



V. CONCLUSIONS 



Jo — 




2|e 2 | 



We have derived the analytic formulas for describing 
the full-time dynamics of the modulational instability in 
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FIG. 2: Full-time evolution of the 4>o component due to modu- 
lational instability. The parameters are the same as in Fig. [T] 



the integrable model of F = 1 Bose-Einstein conden- 
sates. Our results are based on the exact homoclinic 
solution of the matrix NLS equation with the continuous 
plane waves as an initial condition. We have shown that 
there exist cycles of the MI evolution with the reversal 
property when the exponential growth of the modula- 
tion amplitude changes to its exponential decay. The 
solution we present here is an example of large- amplitude 
periodic solutions. It describes an exponential growth of 



a weak modulation of a background for an initial stage 
of the condensate evolution, and in this sense the back- 
ground is unstable. However, in the nonlinear regime this 
exponentially growing mode saturates and subsequently 
transforms into oscillations. As expected, the integrable 
model (|2.3|) does not exhibit long-time chaotic dynam- 
ics contrary to the regimes observed numerically for a 
general case but it may serve as a good analytical 
approximation of the evolution of the condensate experi- 
encing the instability. Higher-order homoclinic solutions 
which correspond to several complex double points can 
be obtained by the method described in [261 ]. 

Strictly speaking, the analysis based on the continu- 
ous wave model is not applicable to the trapped systems. 
Nevertheless, such an approach remains valid when the 
typical spatial extent of the condensate is larger than 
the period of the localized pattern formed in result of 
the instability. More realistic models should account for 
a (small) deviation of the condensate parameters from 
the constraint which provides integrability of the model. 
There exists an approach [3l[ to reveal a persistence 
of the homoclinic orbit when the integrability condition 
breaks, and therefore to establish analytically the exis- 
tence of chaos. This approach is based on the construc- 
tion of the so-called Melnikov function from the quadratic 
products of the Bloch functions evaluated on the homo- 
clinic orbit. In this paper we have explicitly built all 
the ingredients to perform the Melnikov analysis. Corre- 
sponding results will be published elsewhere. 
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